NCC ESCC_Natera analysis 04262024

library(swimplot) library(grid) library(gtable) library(readr) library(mosaic) library(dplyr) library(survival) library(survminer) library(ggplot2) library(scales) library(coxphf) library(ggthemes) library(tidyverse) library(gtsummary) library(flextable) library(parameters) library(car) library(ComplexHeatmap) library(tidyverse) library(readxl) library(survival) library(janitor) library(DT)

#ctDNA positivity by stage and window

#Number of Pts at Baseline - percentage positivity by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.prenac!="",]
circ_datadf <- as.data.frame(circ_data)

total_base <- sum(!is.na(circ_data$ctDNA.prenac))
print(total_base)

[1] 28

circ_data$ctDNA.prenac <- as.factor(circ_data$ctDNA.prenac)
cont_table_base <- table(circ_data$cStage, circ_data$ctDNA.prenac)
print(cont_table_base)

     
      NEGATIVE POSITIVE
  I          1        1
  II         0        3
  III        0       15
  IV         0        8

#Number of Pts post NAT - percentage positivity by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.postnac!="",]
circ_datadf <- as.data.frame(circ_data)

total_postnat <- sum(!is.na(circ_data$ctDNA.postnac))
print(total_postnat)

[1] 27

circ_data$ctDNA.postnac <- as.factor(circ_data$ctDNA.postnac)
cont_table_postnat <- table(circ_data$cStage, circ_data$ctDNA.postnac)
print(cont_table_postnat)

     
      NEGATIVE POSITIVE
  I          2        0
  II         1        2
  III       10        4
  IV         5        3

#Number of Pts post-surgery - percentage positivity by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.4w!="",]
circ_datadf <- as.data.frame(circ_data)

total_mrd <- sum(!is.na(circ_data$ctDNA.4w))
print(total_mrd)

[1] 28

circ_data$ctDNA.4w <- as.factor(circ_data$ctDNA.4w)
cont_table_mrd <- table(circ_data$cStage, circ_data$ctDNA.4w)
print(cont_table_mrd)

     
      NEGATIVE POSITIVE
  I          2        0
  II         2        1
  III       14        1
  IV         7        1

#Summary Table

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_datadf <- as.data.frame(circ_data)

circ_data_subset <- circ_data %>%
  select(
    Gender,
    Age,
    PS,
    cStage,
    TNM.Detailed,
    Neo.Regiment,
    ypStage,
    ypTNM.Detailed,
    TRG,
    pCR,
    Adj.Regimen,
    RFS.Event,
    OS.Event,
    RFS.months,
    OS.months) %>%
  mutate(
    Gender = factor(Gender),
    Age = as.numeric(Age),
    PS = factor(PS),
    cStage = factor(cStage),
    TNM.Detailed = factor(TNM.Detailed),
    Neo.Regiment = factor(Neo.Regiment),
    ypStage = factor(ypStage),
    ypTNM.Detailed = factor(ypTNM.Detailed),
    TRG = factor(TRG),
    pCR = factor(pCR),
    Adj.Regimen = factor(Adj.Regimen),
    RFS.Event = factor(RFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Recurrence", "Recurrence")),
    OS.Event = factor(OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased")),
    RFS.months = as.numeric(RFS.months),
    OS.months = as.numeric(OS.months))
table1 <- circ_data_subset %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
table1

Characteristic	N = 281
Gender
FEMALE	1 (3.6%)
MALE	27 (96%)
Age	71 (50 - 82)
PS
0	26 (93%)
1	2 (7.1%)
cStage
I	2 (7.1%)
II	3 (11%)
III	15 (54%)
IV	8 (29%)
TNM.Detailed
T1-T2N1M0	3 (11%)
T3-T4N2-N3M0	3 (11%)
T3-T4N2-N3M1	4 (14%)
T3N0M0	2 (7.1%)
T3N0M1	1 (3.6%)
T3N1-N2M0	15 (54%)
Neo.Regiment
DCF	19 (68%)
FOLFOX6	6 (21%)
FP	3 (11%)
ypStage
0	4 (14%)
I	3 (11%)
II	8 (29%)
III	7 (25%)
IV	6 (21%)
ypTNM.Detailed
T0-TisN0M0	4 (14%)
T1-T3N0M0	6 (21%)
T1-T3N1-N3M1	5 (18%)
T1-T3N1M0	8 (29%)
T1N0M0	2 (7.1%)
T3N2-N3M0	3 (11%)
TRG
1	12 (43%)
2	13 (46%)
3	3 (11%)
pCR
FALSE	24 (86%)
TRUE	4 (14%)
Adj.Regimen
	25 (89%)
Nivolumab	3 (11%)
RFS.Event
No Recurrence	16 (57%)
Recurrence	12 (43%)
OS.Event
Alive	26 (93%)
Deceased	2 (7.1%)
RFS.months	20 (2 - 36)
OS.months	22 (6 - 34)
1 n (%); Median (Range)

fit1 <- as_flex_table(
  table1,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)

Warning: The `strip_md_bold` argument of `as_flex_table()` is deprecated as of gtsummary 1.6.0.

fit1

Characteristic	N = 281
Gender
FEMALE	1 (3.6%)
MALE	27 (96%)
Age	71 (50 - 82)
PS
0	26 (93%)
1	2 (7.1%)
cStage
I	2 (7.1%)
II	3 (11%)
III	15 (54%)
IV	8 (29%)
TNM.Detailed
T1-T2N1M0	3 (11%)
T3-T4N2-N3M0	3 (11%)
T3-T4N2-N3M1	4 (14%)
T3N0M0	2 (7.1%)
T3N0M1	1 (3.6%)
T3N1-N2M0	15 (54%)
Neo.Regiment
DCF	19 (68%)
FOLFOX6	6 (21%)
FP	3 (11%)
ypStage
0	4 (14%)
I	3 (11%)
II	8 (29%)
III	7 (25%)
IV	6 (21%)
ypTNM.Detailed
T0-TisN0M0	4 (14%)
T1-T3N0M0	6 (21%)
T1-T3N1-N3M1	5 (18%)
T1-T3N1M0	8 (29%)
T1N0M0	2 (7.1%)
T3N2-N3M0	3 (11%)
TRG
1	12 (43%)
2	13 (46%)
3	3 (11%)
pCR
FALSE	24 (86%)
TRUE	4 (14%)
Adj.Regimen
	25 (89%)
Nivolumab	3 (11%)
RFS.Event
No Recurrence	16 (57%)
Recurrence	12 (43%)
OS.Event
Alive	26 (93%)
Deceased	2 (7.1%)
RFS.months	20 (2 - 36)
OS.months	22 (6 - 34)
1n (%); Median (Range)

save_as_docx(fit1, path= "~/Downloads/table1.docx")

#Heatmap for the clinical factors

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_datadf <- as.data.frame(circ_data)
circ_data <- circ_data %>% arrange(cStage)
circ_datadf <- as.data.frame(circ_data)

ha <- HeatmapAnnotation(
  cStage = circ_data$cStage,
  Gender = circ_data$Gender,
  Neo.Regiment = circ_data$Neo.Regiment,
  TRG = circ_data$TRG,
  pCR = circ_data$pCR,
  ypStage = circ_data$ypStage,
  Adj.Regimen = circ_data$Adj.Regimen,
  ctDNA.prenac = circ_data$ctDNA.prenac,
  ctDNA.postnac = circ_data$ctDNA.postnac,
  ctDNA.4w = circ_data$ctDNA.4w,
  ctDNA.12w = circ_data$ctDNA.12w,
  ctDNA.24w = circ_data$ctDNA.24w,
  ctDNA.36w = circ_data$ctDNA.36w,
  ctDNA.48w = circ_data$ctDNA.48w,
  ctDNA.postMRD = circ_data$ctDNA.postMRD,
  RFS.Event = circ_data$RFS.Event,
  OS.Event = circ_data$OS.Event,
  
    col = list(cStage = c("I" = "seagreen1", "II" = "khaki", "III" = "orange", "IV" = "darkmagenta"),
    Gender = c("FEMALE" = "goldenrod" , "MALE" = "blue4"),
    Neo.Regiment = c("FOLFOX6" = "coral", "FP" ="darkgreen", "DCF" ="brown"),
    TRG = c("1" = "yellow3", "2" ="darkgreen", "3" = "brown2"),
    pCR = c("TRUE" = "lightblue", "FALSE" ="orange"),
    ypStage = c("0" = "khaki","I" = "seagreen2", "II" = "cornflowerblue","III" = "orange","IV" ="darkmagenta"),
    Adj.Regimen = c("Nivolumab" = "darkmagenta"),
    ctDNA.prenac = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.postnac = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.4w = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.12w = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.24w = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.36w = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.48w = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    ctDNA.postMRD = c("POSITIVE" = "red3", "NEGATIVE" ="blue"),
    RFS.Event = c("TRUE" = "red3", "FALSE" ="blue"),
    OS.Event = c("TRUE" = "black", "FALSE" ="grey")
)
)
ht <- Heatmap(matrix(nrow = 0, ncol = length(circ_data$cStage)),show_row_names = FALSE,cluster_rows = F,cluster_columns = FALSE, top_annotation = ha)
pdf("heatmap.pdf",width = 7, height = 7)
draw(ht, annotation_legend_side = "bottom")
dev.off()

#Overview plot

setwd("~/Downloads")
clinstage<- read.csv("NCC_ESCC_OP.csv")
clinstage_df<- as.data.frame(clinstage)

#Display the swimmer plot with the label box
oplot<-swimmer_plot(df=clinstage_df,
                    id='PatientName',
                    end='fu.diff.months',
                    fill='gray',
                    width=.01,)
oplot <- oplot + theme(panel.border = element_blank())
oplot <- oplot + scale_y_continuous(breaks = seq(-12, 108, by = 6))
oplot <- oplot + labs(x ="Patients" , y="Months from Surgery")
oplot

##plot events
oplot_ev1 <- oplot + swimmer_points(df_points=clinstage_df,
                                    id='PatientName',
                                    time='date.diff.months',
                                    name_shape ='Event_type',
                                    name_col = 'Event',
                                    size=3.5,fill='black',
                                    #col='darkgreen'
)
oplot_ev1

#Shape customization to Event_type

oplot_ev1.1 <- oplot_ev1 + ggplot2::scale_shape_manual(name="Event_type",values=c(1,16,6,18,18,4),breaks=c('ctDNA_neg','ctDNA_pos', 'Imaging','Surgery','Biopsy', 'Death'))

oplot_ev1.1

#plot treatment

oplot_ev2 <- oplot_ev1.1 + swimmer_lines(df_lines=clinstage_df,
                                         id='PatientName',
                                         start='Tx_start.months',
                                         end='Tx_end.months',
                                         name_col='Tx_type',
                                         size=3.5,
                                         name_alpha = 1.0)
oplot_ev2 <- oplot_ev2 + guides(linetype = guide_legend(override.aes = list(size = 5, color = "black")))
oplot_ev2  

#colour customization
#orange=NAC chemo,orangered=NAC ChemoIO,  purple=ACT, Black=Death, Red=PD, ctDNA negative=white, ctDNA positive=black, cystectomy=blue, TURBT=gray, ACT Chemo=purple, ACT IO=lightblue 
oplot_ev2.2 <- oplot_ev2 + ggplot2::scale_color_manual(name="Event",values=c("lightblue", "gray", "black", "black", "orange", "green", "red", "blue"))
oplot_ev2.2

#RFS by ctDNA post-NAT

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.postnac!="",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.postnac, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.postnac, data = circ_data)

                        n events median 0.95LCL 0.95UCL
ctDNA.postnac=NEGATIVE 18      5     NA      NA      NA
ctDNA.postnac=POSITIVE  9      7   7.85    3.88      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.postnac, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="RFS - ctDNA post-NAC", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.postnac, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.postnac=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     15       3    0.833  0.0878        0.568        0.943
   24      7       2    0.699  0.1141        0.417        0.864

                ctDNA.postnac=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      12.000        4.000        5.000        0.444        0.166        0.136        0.719 

circ_data$ctDNA.postnac <- factor(circ_data$ctDNA.postnac, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.postnac, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.postnac, data = circ_data)

  n= 27, number of events= 12 

                        coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.postnacPOSITIVE 1.5174    4.5603   0.5893 2.575     0.01 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                      exp(coef) exp(-coef) lower .95 upper .95
ctDNA.postnacPOSITIVE      4.56     0.2193     1.437     14.47

Concordance= 0.698  (se = 0.064 )
Likelihood ratio test= 6.58  on 1 df,   p=0.01
Wald test            = 6.63  on 1 df,   p=0.01
Score (logrank) test = 7.95  on 1 df,   p=0.005

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 4.56 (1.44-14.47); p = 0.01"

#RFS by ctDNA Clearance post-NAC - 3 Groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.Clearance!="",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.Clearance, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.Clearance, data = circ_data)

                              n events median 0.95LCL 0.95UCL
ctDNA.Clearance=No Clearance  9      7   7.85    3.88      NA
ctDNA.Clearance=Sustained    13      2     NA      NA      NA
ctDNA.Clearance=Transient     4      3  15.70    8.35      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Clearance, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue","green"), title="RFS - ctDNA Clearance post-NAC", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("No Clearance", "Sustained", "Transient"), legend.title="")

summary(KM_curve, times= c(12, 18, 24))

Call: survfit(formula = surv_object ~ ctDNA.Clearance, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Clearance=No Clearance 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      4       5    0.444   0.166       0.1359        0.719
   18      3       1    0.333   0.157       0.0783        0.623

                ctDNA.Clearance=Sustained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     11       2    0.846     0.1        0.512        0.959
   18      9       0    0.846     0.1        0.512        0.959
   24      7       0    0.846     0.1        0.512        0.959

                ctDNA.Clearance=Transient 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      3       1    0.750   0.217        0.128        0.961
   18      1       1    0.375   0.286        0.011        0.808

circ_data$ctDNA.Clearance <- factor(circ_data$ctDNA.Clearance, levels=c("Sustained","Transient", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Clearance, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Clearance, data = circ_data)

  n= 26, number of events= 12 

                              coef exp(coef) se(coef)     z Pr(>|z|)   
ctDNA.ClearanceTransient    1.9042    6.7139   0.9275 2.053  0.04007 * 
ctDNA.ClearanceNo Clearance 2.1652    8.7160   0.8050 2.690  0.00715 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                            exp(coef) exp(-coef) lower .95 upper .95
ctDNA.ClearanceTransient        6.714     0.1489     1.090     41.35
ctDNA.ClearanceNo Clearance     8.716     0.1147     1.799     42.22

Concordance= 0.747  (se = 0.071 )
Likelihood ratio test= 10.2  on 2 df,   p=0.006
Wald test            = 7.33  on 2 df,   p=0.03
Score (logrank) test = 10.27  on 2 df,   p=0.006

#RFS by ctDNA at 4 weeks

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.4w!="",]

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.4w, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.4w, data = circ_data)

                   n events median 0.95LCL 0.95UCL
ctDNA.4w=NEGATIVE 25      9     NA   19.58      NA
ctDNA.4w=POSITIVE  3      3   2.83    1.87      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.4w, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="RFS - ctDNA at week 4", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

circ_data$ctDNA.4w <- factor(circ_data$ctDNA.4w, levels=c("NEGATIVE","POSITIVE"))
summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.4w, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ctDNA.4w=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     25       0    1.000   0.000        1.000        1.000
   24      8       9    0.618   0.101        0.391        0.781

                ctDNA.4w=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            3            0            1            0            1            1 

cox_fit <- coxph(surv_object ~ ctDNA.4w, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.4w, data = circ_data)

  n= 28, number of events= 12 

                   coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.4wPOSITIVE  2.609    13.585    0.778 3.353 0.000798 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
ctDNA.4wPOSITIVE     13.58    0.07361     2.957     62.42

Concordance= 0.635  (se = 0.062 )
Likelihood ratio test= 8.4  on 1 df,   p=0.004
Wald test            = 11.25  on 1 df,   p=8e-04
Score (logrank) test = 18.68  on 1 df,   p=2e-05

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 13.58 (2.96-62.42); p = 0.001"

#RFS by ctDNA MRD Window (2-16 weeks post-op)

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                    n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 23      7     NA   19.75      NA
ctDNA.MRD=POSITIVE  5      5   3.88    2.83      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="RFS - ctDNA MRD Window", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     23       0    1.000   0.000         1.00        1.000
   24      8       7    0.672   0.103         0.43        0.829

                ctDNA.MRD=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            5            0            1            0            1            1 

cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 28, number of events= 12 

                     coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRDPOSITIVE  3.4335   30.9862   0.8674 3.958 7.54e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     30.99    0.03227      5.66     169.6

Concordance= 0.723  (se = 0.06 )
Likelihood ratio test= 16.8  on 1 df,   p=4e-05
Wald test            = 15.67  on 1 df,   p=8e-05
Score (logrank) test = 32.96  on 1 df,   p=9e-09

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 30.99 (5.66-169.63); p = 0"

#RFS by ctDNA post MRD

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.postMRD!="",]

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.postMRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.postMRD, data = circ_data)

                        n events median 0.95LCL 0.95UCL
ctDNA.postMRD=NEGATIVE 17      3     NA      NA      NA
ctDNA.postMRD=POSITIVE  5      4   8.35    7.85      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.postMRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="RFS - ctDNA Surveillance (post-)MRD", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

circ_data$ctDNA.postMRD <- factor(circ_data$ctDNA.postMRD, levels=c("NEGATIVE","POSITIVE"))
summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.postMRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.postMRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     17       0    1.000  0.0000        1.000        1.000
   24      8       3    0.819  0.0946        0.538        0.938

                ctDNA.postMRD=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            5            0            1            0            1            1 

cox_fit <- coxph(surv_object ~ ctDNA.postMRD, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.postMRD, data = circ_data)

  n= 22, number of events= 7 

                        coef exp(coef) se(coef)     z Pr(>|z|)   
ctDNA.postMRDPOSITIVE  3.177    23.975    1.156 2.748    0.006 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                      exp(coef) exp(-coef) lower .95 upper .95
ctDNA.postMRDPOSITIVE     23.97    0.04171     2.487     231.1

Concordance= 0.75  (se = 0.079 )
Likelihood ratio test= 10.18  on 1 df,   p=0.001
Wald test            = 7.55  on 1 df,   p=0.006
Score (logrank) test = 14.78  on 1 df,   p=1e-04

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 23.97 (2.49-231.13); p = 0.006"

#Time-dependent analysis for post-MRD (>16 weeks) time points

rm(list=ls())
setwd("~/Downloads")
dt <- read_xlsx("NCC ESCC Time Dependent Data_postMRD.xlsx") |>
  clean_names() |>
  mutate(across(.cols = c(window_start_date,dfs_date,
                          surveillance_1_date:surveillance_12_date), 
                .fns = ~ as_date(as.Date(.x, format = "%Y-%m-%d"))))

dt_biomarker <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date,
         surveillance_1_status:surveillance_12_date) |>
  filter(ct_dna_surveillance_available) |>
  pivot_longer(cols = surveillance_1_status:surveillance_12_date,
               names_to = c("visit_number", ".value"),
               names_pattern = "surveillance_(.)_(.*)") |>
  mutate(biomarker_time = day(days(date - window_start_date))) |>
  select(pts_id, biomarker_time, biomarker_status = status) |>
  filter(!is.na(biomarker_time))

glimpse(dt_biomarker)

Rows: 105
Columns: 3
$ pts_id           <chr> "SIG-E003", "SIG-E003", "SIG-E003", "SIG-E003", "SIG-E003", "SIG-E004", "SIG-E004", "SIG-E005", "SIG-E005", "SIG-E00…
$ biomarker_time   <dbl> -90, -40, 50, 141, 222, -102, -44, -90, -27, 57, 141, 218, -98, -49, 119, 245, -91, -34, 50, 144, 222, -100, -30, 54…
$ biomarker_status <chr> "NEGATIVE", "NEGATIVE", "NEGATIVE", "NEGATIVE", "NEGATIVE", "NEGATIVE", "POSITIVE", "NEGATIVE", "NEGATIVE", "NEGATIV…

dt_survival <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date:dfs_date, dfs_event) |>  # Added dfs_event here
  filter(ct_dna_surveillance_available) |>
  mutate(dfs_time = (dfs_date - window_start_date),
         dfs_time = day(days(dfs_time)),
         dfs_event = as.numeric(dfs_event)) |>
  select(pts_id, dfs_time, dfs_event)

glimpse(dt_survival)

Rows: 25
Columns: 3
$ pts_id    <chr> "SIG-E003", "SIG-E004", "SIG-E005", "SIG-E011", "SIG-E013", "SIG-E015", "SIG-E020", "SIG-E021", "SIG-E024", "SIG-E025", "SI…
$ dfs_time  <dbl> 966, -2, 966, 966, 784, 782, 827, 786, 119, -34, 620, 602, 604, -63, 295, 604, 49, 358, 476, 301, 606, 134, 415, 482, 52
$ dfs_event <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1

aux <- dt_survival %>% 
  filter(dfs_time <= 0)

tab <- left_join(aux, dt) |>
  select(pts_id, window_start_date, dfs_time, dfs_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = dfs_date:surveillance_12_date, 
                .fns = ~ as_date(.x))) |>
  select(pts_id, window_start_date, dfs_date, dfs_time)

Joining with `by = join_by(pts_id, dfs_event)`

datatable(tab, filter = "top")

dt_survival <- dt_survival |>
  filter(dfs_time > 0)

aux <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date, dfs_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = surveillance_1_date:surveillance_12_date, 
                .fns = ~ .x - window_start_date)) |>
  mutate(across(.cols = surveillance_1_date:surveillance_12_date, 
                .fns = ~ .x < 0)) |>
  rowwise() |>
  mutate(sum_neg = 
           sum(c_across(surveillance_1_date:surveillance_12_date),
               na.rm = TRUE))  |>
  select(pts_id, sum_neg)

tab <- left_join(aux, dt) |>
  filter(sum_neg > 0) |>
  select(pts_id, sum_neg, window_start_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = window_start_date:surveillance_12_date, 
                .fns = ~ as_date(.x))) 

Joining with `by = join_by(pts_id)`

datatable(tab, filter = "top")

aux <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date, dfs_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = dfs_date:surveillance_12_date, 
                .fns = ~ .x - window_start_date)) |>
  mutate(across(.cols = surveillance_2_date:surveillance_12_date,
                .fns = ~ dfs_date < .x)) |>
  rowwise() |>
  mutate(n_biomarker_after_event = sum(c_across(surveillance_2_date:
                                                  surveillance_12_date), 
                                       na.rm = TRUE)) |>
  mutate(across(.cols = surveillance_1_date:surveillance_12_date,
                .fns = ~ !is.na(.x))) |>
  mutate(total_biomarker = sum(c_across(surveillance_2_date:
                                          surveillance_12_date), 
                               na.rm = TRUE)) |>
  select(pts_id, n_biomarker_after_event, total_biomarker)

temp <- aux |> 
  select(-pts_id) |>
  group_by(n_biomarker_after_event, total_biomarker) |>
  summarize(freq = n())

`summarise()` has grouped output by 'n_biomarker_after_event'. You can override using the `.groups` argument.

tab <- left_join(aux, dt) |>
  select(pts_id, n_biomarker_after_event, total_biomarker, 
         dfs_date,
         surveillance_2_date:surveillance_12_date) |>
  mutate(across(.cols = dfs_date:surveillance_12_date, 
                .fns = ~ as_date(.x))) |>
  filter(n_biomarker_after_event > 0)

Joining with `by = join_by(pts_id)`

datatable(tab, filter = "top")

aux <- tmerge(data1 = dt_survival, 
              data2 = dt_survival,
              id = pts_id, 
              dfs_event = event(dfs_time, dfs_event))
dt_final <- tmerge(data1 = aux, 
                   data2 = dt_biomarker,
                   id = pts_id, 
                   biomarker_status = 
                     tdc(biomarker_time, biomarker_status))

datatable(dt_final, filter = "top")

# Syntax if there is not time-dependent covariate
# fit <- coxph(Surv(dfs_time, dfs_event) ~ biomarker_status,
#              data = dt_final)
# summary(fit)

fit <- coxph(Surv(tstart, tstop, dfs_event) ~ biomarker_status,
             data = dt_final)
summary(fit)

Call:
coxph(formula = Surv(tstart, tstop, dfs_event) ~ biomarker_status, 
    data = dt_final)

  n= 76, number of events= 7 

                           coef exp(coef) se(coef)     z Pr(>|z|)   
biomarker_statusPOSITIVE  3.433    30.965    1.102 3.115  0.00184 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                         exp(coef) exp(-coef) lower .95 upper .95
biomarker_statusPOSITIVE     30.96    0.03229     3.571     268.5

Concordance= 0.84  (se = 0.083 )
Likelihood ratio test= 15.55  on 1 df,   p=8e-05
Wald test            = 9.7  on 1 df,   p=0.002
Score (logrank) test = 21.42  on 1 df,   p=4e-06

#RFS by ctDNA anytime post-surgery

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.anytime!="",]

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.anytime, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.anytime, data = circ_data)

                        n events median 0.95LCL 0.95UCL
ctDNA.anytime=NEGATIVE 17      2     NA      NA      NA
ctDNA.anytime=POSITIVE 11     10   8.35    5.55      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.anytime, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="RFS - ctDNA anytime post-surgery", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

circ_data$ctDNA.anytime <- factor(circ_data$ctDNA.anytime, levels=c("NEGATIVE","POSITIVE"))
summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.anytime, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.anytime=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     17       0    1.000  0.0000        1.000        1.000
   24      8       2    0.882  0.0781        0.606        0.969

                ctDNA.anytime=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           11            0            1            0            1            1 

cox_fit <- coxph(surv_object ~ ctDNA.anytime, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.anytime, data = circ_data)

  n= 28, number of events= 12 

                         coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.anytimePOSITIVE  2.8517   17.3168   0.8012 3.559 0.000372 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                      exp(coef) exp(-coef) lower .95 upper .95
ctDNA.anytimePOSITIVE     17.32    0.05775     3.601     83.27

Concordance= 0.787  (se = 0.056 )
Likelihood ratio test= 18.5  on 1 df,   p=2e-05
Wald test            = 12.67  on 1 df,   p=4e-04
Score (logrank) test = 21.4  on 1 df,   p=4e-06

cox_fit_summary <- summary(cox_fit)

#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 17.32 (3.6-83.27); p = 0"

#Time-dependent analysis for post-surgery time points

rm(list=ls())
setwd("~/Downloads")
dt <- read_xlsx("NCC ESCC Time Dependent Data_anytime postop.xlsx") |>
  clean_names() |>
  mutate(across(.cols = c(window_start_date,dfs_date,
                          surveillance_1_date:surveillance_12_date), 
                .fns = ~ as_date(as.Date(.x, format = "%Y-%m-%d"))))

dt_biomarker <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date,
         surveillance_1_status:surveillance_12_date) |>
  filter(ct_dna_surveillance_available) |>
  pivot_longer(cols = surveillance_1_status:surveillance_12_date,
               names_to = c("visit_number", ".value"),
               names_pattern = "surveillance_(.)_(.*)") |>
  mutate(biomarker_time = day(days(date - window_start_date))) |>
  select(pts_id, biomarker_time, biomarker_status = status) |>
  filter(!is.na(biomarker_time))

glimpse(dt_biomarker)

Rows: 117
Columns: 3
$ pts_id           <chr> "SIG-E003", "SIG-E003", "SIG-E003", "SIG-E003", "SIG-E003", "SIG-E004", "SIG-E004", "SIG-E005", "SIG-E005", "SIG-E00…
$ biomarker_time   <dbl> 30, 80, 170, 261, 342, 18, 76, 30, 93, 177, 261, 338, 22, 71, 239, 365, 29, 86, 170, 264, 342, 20, 90, 174, 258, 349…
$ biomarker_status <chr> "NEGATIVE", "NEGATIVE", "NEGATIVE", "NEGATIVE", "NEGATIVE", "NEGATIVE", "POSITIVE", "NEGATIVE", "NEGATIVE", "NEGATIV…

dt_survival <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date:dfs_date, dfs_event) |>  # Added dfs_event here
  filter(ct_dna_surveillance_available) |>
  mutate(dfs_time = (dfs_date - window_start_date),
         dfs_time = day(days(dfs_time)),
         dfs_event = as.numeric(dfs_event)) |>
  select(pts_id, dfs_time, dfs_event)

glimpse(dt_survival)

Rows: 28
Columns: 3
$ pts_id    <chr> "SIG-E003", "SIG-E004", "SIG-E005", "SIG-E011", "SIG-E013", "SIG-E015", "SIG-E020", "SIG-E021", "SIG-E024", "SIG-E025", "SI…
$ dfs_time  <dbl> 1086, 118, 1086, 1086, 904, 902, 947, 906, 239, 86, 740, 722, 724, 57, 415, 724, 169, 478, 596, 421, 726, 254, 535, 602, 34…
$ dfs_event <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1

aux <- dt_survival %>% 
  filter(dfs_time <= 0)

tab <- left_join(aux, dt) |>
  select(pts_id, window_start_date, dfs_time, dfs_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = dfs_date:surveillance_12_date, 
                .fns = ~ as_date(.x))) |>
  select(pts_id, window_start_date, dfs_date, dfs_time)

Joining with `by = join_by(pts_id, dfs_event)`

datatable(tab, filter = "top")

dt_survival <- dt_survival |>
  filter(dfs_time > 0)

aux <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date, dfs_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = surveillance_1_date:surveillance_12_date, 
                .fns = ~ .x - window_start_date)) |>
  mutate(across(.cols = surveillance_1_date:surveillance_12_date, 
                .fns = ~ .x < 0)) |>
  rowwise() |>
  mutate(sum_neg = 
           sum(c_across(surveillance_1_date:surveillance_12_date),
               na.rm = TRUE))  |>
  select(pts_id, sum_neg)

tab <- left_join(aux, dt) |>
  filter(sum_neg > 0) |>
  select(pts_id, sum_neg, window_start_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = window_start_date:surveillance_12_date, 
                .fns = ~ as_date(.x))) 

Joining with `by = join_by(pts_id)`

datatable(tab, filter = "top")

aux <- dt |>
  select(pts_id, ct_dna_surveillance_available,
         window_start_date, dfs_date,
         surveillance_1_date:surveillance_12_date) |>
  mutate(across(.cols = dfs_date:surveillance_12_date, 
                .fns = ~ .x - window_start_date)) |>
  mutate(across(.cols = surveillance_2_date:surveillance_12_date,
                .fns = ~ dfs_date < .x)) |>
  rowwise() |>
  mutate(n_biomarker_after_event = sum(c_across(surveillance_2_date:
                                                  surveillance_12_date), 
                                       na.rm = TRUE)) |>
  mutate(across(.cols = surveillance_1_date:surveillance_12_date,
                .fns = ~ !is.na(.x))) |>
  mutate(total_biomarker = sum(c_across(surveillance_2_date:
                                          surveillance_12_date), 
                               na.rm = TRUE)) |>
  select(pts_id, n_biomarker_after_event, total_biomarker)

temp <- aux |> 
  select(-pts_id) |>
  group_by(n_biomarker_after_event, total_biomarker) |>
  summarize(freq = n())

`summarise()` has grouped output by 'n_biomarker_after_event'. You can override using the `.groups` argument.

tab <- left_join(aux, dt) |>
  select(pts_id, n_biomarker_after_event, total_biomarker, 
         dfs_date,
         surveillance_2_date:surveillance_12_date) |>
  mutate(across(.cols = dfs_date:surveillance_12_date, 
                .fns = ~ as_date(.x))) |>
  filter(n_biomarker_after_event > 0)

Joining with `by = join_by(pts_id)`

datatable(tab, filter = "top")

aux <- tmerge(data1 = dt_survival, 
              data2 = dt_survival,
              id = pts_id, 
              dfs_event = event(dfs_time, dfs_event))
dt_final <- tmerge(data1 = aux, 
                   data2 = dt_biomarker,
                   id = pts_id, 
                   biomarker_status = 
                     tdc(biomarker_time, biomarker_status))

datatable(dt_final, filter = "top")

# Syntax if there is not time-dependent covariate
# fit <- coxph(Surv(dfs_time, dfs_event) ~ biomarker_status,
#              data = dt_final)
# summary(fit)

fit <- coxph(Surv(tstart, tstop, dfs_event) ~ biomarker_status,
             data = dt_final)
summary(fit)

Call:
coxph(formula = Surv(tstart, tstop, dfs_event) ~ biomarker_status, 
    data = dt_final)

  n= 113, number of events= 12 
   (28 observations deleted due to missingness)

                            coef exp(coef) se(coef)     z Pr(>|z|)    
biomarker_statusPOSITIVE  3.2855   26.7230   0.7845 4.188 2.81e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                         exp(coef) exp(-coef) lower .95 upper .95
biomarker_statusPOSITIVE     26.72    0.03742     5.743     124.3

Concordance= 0.843  (se = 0.059 )
Likelihood ratio test= 26.03  on 1 df,   p=3e-07
Wald test            = 17.54  on 1 df,   p=3e-05
Score (logrank) test = 38  on 1 df,   p=7e-10

#RFS by ctDNA Dynamics post-NAC to Week 4 - 3 groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.postnac == "NEGATIVE" & ctDNA.4w == "NEGATIVE" ~ "1",
    ctDNA.postnac == "POSITIVE" & ctDNA.4w == "NEGATIVE" ~ "2",
    ctDNA.postnac == "POSITIVE" & ctDNA.4w == "POSITIVE" ~ "3"
  ))

circ_data <- circ_data[!is.na(circ_data$ctDNA.Dynamics),]
survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 17      4     NA      NA      NA
ctDNA.Dynamics=2  7      5  13.83    5.55      NA
ctDNA.Dynamics=3  2      2   2.35    1.87      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","red"), title="RFS - ctDNA Dynamics post-NAC - Week 4", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("Persistently Negative","Converted Negative", "Persistently Positive"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     17       0     1.00   0.000        1.000        1.000
   24      7       4     0.74   0.113        0.443        0.895

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            7            0            1            0            1            1 

                ctDNA.Dynamics=3 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            2            0            1            0            1            1 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2","3"), labels=c("Persistently Negative","Converted Negative", "Persistently Positive"))
cox_fit <- coxphf(surv_object ~ ctDNA.Dynamics, data=circ_data)
summary(cox_fit)

coxphf(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                                        coef se(coef)  exp(coef) lower 0.95  upper 0.95    Chisq          p
ctDNA.DynamicsConverted Negative    1.415181 0.660610   4.117232   1.162499    15.37828  4.76164 0.02910081
ctDNA.DynamicsPersistently Positive 5.119773 1.952631 167.297395  11.566163 24434.33628 13.90506 0.00019228

Likelihood ratio test=16.3327 on 2 df, p=0.000284053, n=26
Wald test = 9.50381 on 2 df, p = 0.00863523

Covariance-Matrix:
                                    ctDNA.DynamicsConverted Negative ctDNA.DynamicsPersistently Positive
ctDNA.DynamicsConverted Negative                           0.4364055                           0.2744938
ctDNA.DynamicsPersistently Positive                        0.2744938                           3.8127677

#RFS by ctDNA Dynamics post-NAC to Week 4 - 4 groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.postnac == "NEGATIVE" & ctDNA.4w == "NEGATIVE" ~ "1",
    ctDNA.postnac == "POSITIVE" & ctDNA.4w == "NEGATIVE" ~ "2",
    ctDNA.postnac == "NEGATIVE" & ctDNA.4w == "POSITIVE" ~ "3",
    ctDNA.postnac == "POSITIVE" & ctDNA.4w == "POSITIVE" ~ "4"
  ))

circ_data <- circ_data[!is.na(circ_data$ctDNA.Dynamics),]
survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 17      4     NA      NA      NA
ctDNA.Dynamics=2  7      5  13.83    5.55      NA
ctDNA.Dynamics=3  1      1   8.35      NA      NA
ctDNA.Dynamics=4  2      2   2.35    1.87      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple","red"), title="RFS - ctDNA Dynamics post-NAC - Week 4", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("Persistently Negative","Converted Negative", "Converted Positive","Persistently Positive"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     17       0     1.00   0.000        1.000        1.000
   24      7       4     0.74   0.113        0.443        0.895

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            7            0            1            0            1            1 

                ctDNA.Dynamics=3 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            1            0            1            0            1            1 

                ctDNA.Dynamics=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            2            0            1            0            1            1 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2","3","4"), labels=c("Persistently Negative","Converted Negative","Converted Positive","Persistently Positive"))
cox_fit <- coxphf(surv_object ~ ctDNA.Dynamics, data=circ_data)
summary(cox_fit)

coxphf(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                                        coef  se(coef)  exp(coef) lower 0.95  upper 0.95     Chisq            p
ctDNA.DynamicsConverted Negative    1.440334 0.6697529   4.222106   1.192923    15.76026  4.928230 0.0264213641
ctDNA.DynamicsConverted Positive    2.517717 1.0853732  12.400250   1.153475    81.17250  4.214575 0.0400781097
ctDNA.DynamicsPersistently Positive 5.319678 1.9341003 204.318068  13.851382 30074.62703 14.722641 0.0001245417

Likelihood ratio test=19.40201 on 3 df, p=0.0002257542, n=27
Wald test = 12.0956 on 3 df, p = 0.007062756

Covariance-Matrix:
                                    ctDNA.DynamicsConverted Negative ctDNA.DynamicsConverted Positive ctDNA.DynamicsPersistently Positive
ctDNA.DynamicsConverted Negative                           0.4485690                        0.2580369                           0.2791462
ctDNA.DynamicsConverted Positive                           0.2580369                        1.1780349                           0.3771376
ctDNA.DynamicsPersistently Positive                        0.2791462                        0.3771376                           3.7407440

#RFS by ctDNA Dynamics post-NAC to MRD - 4 groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.postnac == "NEGATIVE" & ctDNA.MRD == "NEGATIVE" ~ "1",
    ctDNA.postnac == "POSITIVE" & ctDNA.MRD == "NEGATIVE" ~ "2",
    ctDNA.postnac == "NEGATIVE" & ctDNA.MRD == "POSITIVE" ~ "3",
    ctDNA.postnac == "POSITIVE" & ctDNA.MRD == "POSITIVE" ~ "4"
  ))

circ_data <- circ_data[!is.na(circ_data$ctDNA.Dynamics),]
survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 17      4     NA      NA      NA
ctDNA.Dynamics=2  5      3  19.75   13.83      NA
ctDNA.Dynamics=3  1      1   8.35      NA      NA
ctDNA.Dynamics=4  4      4   3.35    1.87      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple","red"), title="RFS - ctDNA Dynamics post-NAC - MRD Window", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", legend.labs=c("Persistently Negative","Converted Negative", "Converted Positive","Persistently Positive"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     17       0     1.00   0.000        1.000        1.000
   24      7       4     0.74   0.113        0.443        0.895

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            5            0            1            0            1            1 

                ctDNA.Dynamics=3 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            1            0            1            0            1            1 

                ctDNA.Dynamics=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            4            0            1            0            1            1 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2","3","4"), labels=c("Persistently Negative","Converted Negative","Converted Positive","Persistently Positive"))
cox_fit <- coxphf(surv_object ~ ctDNA.Dynamics, data=circ_data)
summary(cox_fit)

coxphf(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                                        coef  se(coef)  exp(coef) lower 0.95  upper 0.95     Chisq            p
ctDNA.DynamicsConverted Negative    1.050415 0.7518683   2.858836  0.6387374    11.75356  1.993070 1.580203e-01
ctDNA.DynamicsConverted Positive    2.819645 1.1734372  16.770898  1.4543670   140.18499  4.827784 2.800461e-02
ctDNA.DynamicsPersistently Positive 5.277803 1.7038744 195.939019 16.5092118 28090.67912 21.337835 3.850563e-06

Likelihood ratio test=23.80561 on 3 df, p=2.742645e-05, n=27
Wald test = 12.33767 on 3 df, p = 0.006311568

Covariance-Matrix:
                                    ctDNA.DynamicsConverted Negative ctDNA.DynamicsConverted Positive ctDNA.DynamicsPersistently Positive
ctDNA.DynamicsConverted Negative                           0.5653059                        0.2539661                           0.2567316
ctDNA.DynamicsConverted Positive                           0.2539661                        1.3769549                           0.5559755
ctDNA.DynamicsPersistently Positive                        0.2567316                        0.5559755                           2.9031880

#Association of ctDNA post-NAC and TRG status

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")

#Vertical Fisher plot for ctDNA clearance post-NAT and Rec Status
circ_data <- circ_data[!is.na(circ_data$ctDNA.postnac),]
circ_data$ctDNA.postnac <- factor(circ_data$ctDNA.postnac, levels=c("NEGATIVE","POSITIVE"))
circ_data$TRG <- factor(circ_data$TRG, levels = c("1", "2", "3"), labels = c("TRG1", "TRG2", "TRG3"))
contingency_table <- table(circ_data$TRG, circ_data$ctDNA.postnac)
fisher_exact_test <- fisher.test(contingency_table)
chi_square_test <- chisq.test(contingency_table)

Warning: Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 13.028, df = 2, p-value = 0.001483

print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.001404
alternative hypothesis: two.sided

print(contingency_table)

      
       NEGATIVE POSITIVE
  TRG1        3        8
  TRG2       12        1
  TRG3        3        0

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2

# Swapping x and y in ggplot function to make bar plot vertical
ggplot(table_df, aes(y = Var1, x = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(x = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "ctDNA post-NAT & TRG", y = "TRG Status", x = "Patients (%)", fill = "ctDNA post-NAC") +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("NEGATIVE" = "lightblue3", "POSITIVE" = "red")) + # define custom colors
  theme(axis.text.y = element_text(angle = 0, hjust = 1.5, size = 14), # increase y-axis text size
        axis.text.x = element_text(size = 14, color = "black"), # increase x-axis text size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#RFS by ctDNA post-NAC and TRG - 3 groups

setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.postnac!="",]

circ_data$ctDNA.pCR <- NA
circ_data <- circ_data %>%
  mutate(ctDNA.pCR = case_when(
    ctDNA.postnac == "NEGATIVE" & TRG %in% c(1, 2) ~ "1",
    ctDNA.postnac == "NEGATIVE" & TRG == 3 ~ "2",
    ctDNA.postnac == "POSITIVE" & TRG %in% c(1, 2) ~ "3"
  ))

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.pCR, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.pCR, data = circ_data)

             n events median 0.95LCL 0.95UCL
ctDNA.pCR=1 15      4     NA      NA      NA
ctDNA.pCR=2  3      1     NA    5.65      NA
ctDNA.pCR=3  9      7   7.85    3.88      NA

circ_data <- circ_data[!is.na(circ_data$ctDNA.pCR),]
circ_data$ctDNA.pCR <- factor(circ_data$ctDNA.pCR, levels=c("1","2", "3"), labels=c("TRG1/2 ctDNA (-)","TRG3 ctDNA (-)", "TRG1/2 ctDNA (+)"))
survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.pCR, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.pCR, data = circ_data)

                            n events median 0.95LCL 0.95UCL
ctDNA.pCR=TRG1/2 ctDNA (-) 15      4     NA      NA      NA
ctDNA.pCR=TRG3 ctDNA (-)    3      1     NA    5.65      NA
ctDNA.pCR=TRG1/2 ctDNA (+)  9      7   7.85    3.88      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.pCR, data = circ_data,conf.int=0.95,conf.type="log-log")
summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.pCR, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.pCR=TRG1/2 ctDNA (-) 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     15       0    1.000   0.000        1.000        1.000
   24      7       4    0.709   0.124        0.395        0.881

                ctDNA.pCR=TRG3 ctDNA (-) 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            3            0            1            0            1            1 

                ctDNA.pCR=TRG1/2 ctDNA (+) 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            9            0            1            0            1            1 

ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, 
           break.time.by=6, palette=c("blue", "green", "red"), title="DFS - ctDNA post-NAC & TRG", ylab= "Recurrence-Free Survival", xlab="Months from Surgery", 
           legend.labs=c("TRG1/2 ctDNA (-)","TRG3 ctDNA (-)", "TRG1/2 ctDNA (+)"), legend.title="")

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.postnac!="",]

circ_data$ctDNA.pCR <- NA
circ_data <- circ_data %>%
  mutate(ctDNA.pCR = case_when(
    ctDNA.postnac == "NEGATIVE" & TRG %in% c(1, 2) ~ "1",
    ctDNA.postnac == "NEGATIVE" & TRG == 3 ~ "2",
    ctDNA.postnac == "POSITIVE" & TRG %in% c(1, 2) ~ "3"
  ))
circ_data$ctDNA.pCR <- factor(circ_data$ctDNA.pCR, levels=c("1","2","3"))
surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.pCR, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.pCR, data = circ_data)

  n= 27, number of events= 12 

             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.pCR2 0.4458    1.5617   1.1200 0.398   0.6906  
ctDNA.pCR3 1.5927    4.9168   0.6308 2.525   0.0116 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.pCR2     1.562     0.6403    0.1739     14.03
ctDNA.pCR3     4.917     0.2034    1.4281     16.93

Concordance= 0.715  (se = 0.064 )
Likelihood ratio test= 6.72  on 2 df,   p=0.03
Wald test            = 6.67  on 2 df,   p=0.04
Score (logrank) test = 8.04  on 2 df,   p=0.02

#RFS by ctDNA post-NAC and pCR - 3 groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.4w!="",]
circ_data$RFS.months=circ_data$RFS.months-1
circ_data <- circ_data[circ_data$RFS.months>=0,]

circ_data$ctDNA.pCR <- NA
circ_data <- circ_data %>%
  mutate(ctDNA.pCR = case_when(
    ctDNA.postnac == "NEGATIVE" & pCR == "TRUE" ~ "1",
    ctDNA.postnac == "NEGATIVE" & pCR == "FALSE" ~ "2",
    ctDNA.postnac == "POSITIVE" & pCR == "FALSE" ~ "3",
  ))

survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.pCR, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.pCR, data = circ_data)

   1 observation deleted due to missingness 
             n events median 0.95LCL 0.95UCL
ctDNA.pCR=1  4      1     NA    4.65      NA
ctDNA.pCR=2 14      4     NA   18.58      NA
ctDNA.pCR=3  9      7   6.85    2.88      NA

circ_data <- circ_data[!is.na(circ_data$ctDNA.pCR),]
circ_data$ctDNA.pCR <- factor(circ_data$ctDNA.pCR, levels=c("1","2", "3"), labels=c("pCR ctDNA (-)","No pCR ctDNA (-)", "No pCR ctDNA (+)"))
survfit(Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)~ctDNA.pCR, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) ~ 
    ctDNA.pCR, data = circ_data)

                            n events median 0.95LCL 0.95UCL
ctDNA.pCR=pCR ctDNA (-)     4      1     NA    4.65      NA
ctDNA.pCR=No pCR ctDNA (-) 14      4     NA   18.58      NA
ctDNA.pCR=No pCR ctDNA (+)  9      7   6.85    2.88      NA

surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.pCR, data = circ_data,conf.int=0.95,conf.type="log-log")
summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.pCR, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.pCR=pCR ctDNA (-) 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0      4       0     1.00   0.000        1.000        1.000
   24      1       1     0.75   0.217        0.128        0.961

                ctDNA.pCR=No pCR ctDNA (-) 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     14       0    1.000   0.000        1.000         1.00
   24      5       4    0.686   0.132        0.359         0.87

                ctDNA.pCR=No pCR ctDNA (+) 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            9            0            1            0            1            1 

ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, 
           break.time.by=6, palette=c("blue", "green", "red"), title="DFS - ctDNA post-NAC & pCR", ylab= "Recurrence-Free Survival", xlab="Time (Months)", 
           legend.labs=c("pCR ctDNA (-)","No pCR ctDNA (-)", "No pCR ctDNA (+)"), legend.title="")

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.4w!="",]
circ_data$RFS.months=circ_data$RFS.months-1
circ_data <- circ_data[circ_data$RFS.months>=0,]

circ_data$ctDNA.pCR <- NA
circ_data <- circ_data %>%
  mutate(ctDNA.pCR = case_when(
    ctDNA.postnac == "NEGATIVE" & pCR == "TRUE" ~ "1",
    ctDNA.postnac == "NEGATIVE" & pCR == "FALSE" ~ "2",
    ctDNA.postnac == "POSITIVE" & pCR == "FALSE" ~ "3",
  ))
circ_data$ctDNA.pCR <- factor(circ_data$ctDNA.pCR, levels=c("1","2","3"))
surv_object <-Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.pCR, data=circ_data)
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.pCR, data = circ_data)

  n= 27, number of events= 12 
   (1 observation deleted due to missingness)

              coef exp(coef) se(coef)     z Pr(>|z|)
ctDNA.pCR2 0.03296   1.03351  1.11864 0.029    0.976
ctDNA.pCR3 1.54364   4.68162  1.07082 1.442    0.149

           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.pCR2     1.034     0.9676    0.1154     9.258
ctDNA.pCR3     4.682     0.2136    0.5740    38.184

Concordance= 0.692  (se = 0.078 )
Likelihood ratio test= 6.58  on 2 df,   p=0.04
Wald test            = 6.63  on 2 df,   p=0.04
Score (logrank) test = 7.95  on 2 df,   p=0.02

#Multivariate cox regression for RFS - ctDNA post-NAC

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.postnac!="",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.postnac <- factor(circ_data$ctDNA.postnac, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels=c("1","2"), labels = c("≤70", ">70"))
circ_data$Stage <- factor(circ_data$Stage, levels = c("I-II", "III-IV"))
surv_object <- Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.postnac + Age.Group + Stage, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for RFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression for RFS - ctDNA week 4

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.4w!="",]

circ_data$ctDNA.4w <- factor(circ_data$ctDNA.4w, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels=c("1","2"), labels = c("≤70", ">70"))
circ_data$Stage <- factor(circ_data$Stage, levels = c("I-II", "III-IV"))
surv_object <- Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.4w + Age.Group + Stage, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for RFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression for RFS - ctDNA MRD Window

rm(list=ls())
setwd("~/Downloads") 
circ_data <- read.csv("NCC ESCC Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels=c("1","2"), labels = c("≤70", ">70"))
circ_data$Stage <- factor(circ_data$Stage, levels = c("I-II", "III-IV"))
surv_object <- Surv(time = circ_data$RFS.months, event = circ_data$RFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.MRD + Age.Group + Stage, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for RFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)